#include "ThresholdEntropy.h"
#define IMG_PXL(__image__, __offset__) ((uchar*)(__image__)->imageData)[__offset__]


ThresholdEntropy::ThresholdEntropy ()
{
	int histSize = 256;
	_histogram = cvCreateHist (1, &histSize, CV_HIST_ARRAY);
}


ThresholdEntropy::~ThresholdEntropy()
{
	cvReleaseHist (&_histogram);
}

void ThresholdEntropy::process (IplImage* inImage, IplImage* outImage, const IplImage* maskImage)
{
	int threshValue = calcThreshold (inImage, maskImage);
	for (int v = 0, lineOffset = 0; v < inImage->height; ++v, lineOffset += inImage->widthStep) {
		for (int u = 0; u < inImage->width; ++u) {
			if (maskImage == NULL || IMG_PXL (maskImage, lineOffset + u) != 0) {
				if (IMG_PXL (inImage, lineOffset + u) > threshValue)
					IMG_PXL (outImage, lineOffset + u) = 255;
				else
					IMG_PXL (outImage, lineOffset + u) = 0;
			}
		}
	}
}

int ThresholdEntropy::calcThreshold (IplImage* inImage,  const IplImage* maskImage)
{
#define HISTOGRAM_BIN(__index__) cvGetReal1D (_histogram->bins, __index__)

	cvCalcHist (&inImage, _histogram, 0, maskImage);

	double gMax = 0;
	int tMax = 0;
	double c0 = 0;

	double c1 = 0;
	for (int i = 0; i < 256; ++i)
		c1 += HISTOGRAM_BIN (i);

	for (int t = 0; t < 256; ++t) {
		double g = 0;
		if (c0 != 0) {
			for (int i = 0; i < t - 1; ++i) {
				if (cvGetReal1D (_histogram->bins, i) != 0) {
					double p = (double)HISTOGRAM_BIN (i) / c0;
					g -= p * log (p);
				}
			}
		}
		if (c1 != 0) {
			for (int i = t; i < 256; ++i) {
				if (HISTOGRAM_BIN(i) != 0) {
					double p = (double)HISTOGRAM_BIN(i) / c1;
					g -= p * log (p);
				}
			}
		}
		c0 += HISTOGRAM_BIN(t);
		c1 -= HISTOGRAM_BIN(t);
		if (t == 0) {
			gMax = g;
			tMax = 0;
		} else if (g > gMax) {
			gMax = g;
			tMax = t;
		}
	}
	return tMax;
}
